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Abstract 

We introduce a new method for reconstructing the primordial power spectrum, P(k), directly 
from observations of the Cosmic Microwave Background (CMB). We employ Singular Value De- 
composition (SVD) to invert the radiation perturbation transfer function. The degeneracy of the 
multipole i to wavenumber k linear mapping is thus reduced. This enables the inversion to be car- 
ried out at each point along a Monte Carlo Markov Chain (MCMC) exploration of the combined 
P(k) and cosmological parameter space. We present best-fit P{k) obtained with this method along 
with other cosmological parameters. 
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I. INTRODUCTION 



— * 

The primordial power spectrum of scalar, curvature perturbations <&(k) is defined as, 

P(k) = ^6\k-k>)mk)$*(k')), (1) 

— * 

where k = \k\ is the wavenumber. The spectrum encodes the initial conditions for the 
system of coupled Einstein-Boltzmann equations which describe the evolution of density 
and radiation perturbations about the FRW background. The spectrum itself is considered 
as a unique window into the era approaching the Planck time. In most models a period 
of scalar field driven inflation is used to solve the cosmological problems and set the near 
scale-invariant form of P(k) prior to the radiation dominated epoch. In this case a simple 
power law parametrised by an amplitude A s and spectral index n s suffices to describe the 
initial conditions to sufficient accuracy for the current data. Higher order contributions such 
as a mild curvature dn s /d\ak have also been explored although it is not strictly required 
by the best-fit models. 

A separate approach to studying the physics behind this early phase is to drop any 
model dependent assumptions of near scale invariance and allow more general functional 
forms of the initial spectrum. One drawback of this approach is the increase in parameter 
space which needs to be explored which increases the complexity of the data fitting step. 
A second drawback is the limited information content of the observations; the effect of 
sample variance, limited range of scales probed, and the degeneracy in the mapping of k 
to angular multipoles on the sky means that one cannot expect to constrain arbitrarily 
complex functions with many degrees of freedom. However there is scope to go beyond the 
model dependent approach. This is particularly true if one is interested in constraining the 
presence of features on the spectrum. These could be in the form of 'glitches' or step-like 
features which would be otherwise unconstrained. In fact many models have been proposed 
which predict features on the spectrum due to, for exam ple, features on the inflaton potential 
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There are two approaches to reconstructing P(k); parameterisation and direct inversion. 
None of the various methods have shown conclusive evidence for a departure from near scale- 
invariance of P{k). Despite this there have been tantalising hints of anomalous features in 
the CMB. One example of this is that the first year WMAP results gave an indication of a 



cut-off in P(k) on large scales. With subsequent data releases the significance of this feature 
has been reduced, although future observations of the polarisation of the CMB may provide 
more conclusive evidence 
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171 ]. In [18j we also showed evidence for a dip in power at 



k » 0.002MPC- 1 . 

Numerous parametertic searches for features with a similar form to those in complex 
inflationary models have been performed along with simple binning of P(k) 
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28j |. Methods of direct inversion which make no assumptions about 



the early universe model being probed [18|, 
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43j are hampered by the singular nature of the transfer function that takes 



P(k) and transfers it onto the CMB or LSS modes. In general this causes the process of 
estimation to be prohibitively slow so as not to allow joint estimation of a free P(k) with 
cosmological parameters. Instead one usually assumes a set of cosmological parameters, 
this allows the use of a fiducial CMB photon transfer function to integrate the primordial 
curvature perturbation into today's photon distribution perturbation. However, this hides a 
significant degeneracy between features in the primordial power spectrum, and the physical 
parameters which determine the height and position of acoustic peaks in the CMB. It has 



been pointed out 



451 ] that adding polarisation information or LSS data can help 



break some of the degeneracies. 

In this paper we propose a new method for direct inversion of P(k) from CMB observation 
using Singular Value Decomposition (SVD). The method is fast enough to allow us to carry 
out a joint estimation of P{k) and the cosmological parameters. The form of P{k) is derived 
from the SVD inversion. We also show that there are regions of k for which polarisation 
data has the potential to more accurately constrain P(k). The paper is organised as follows; 
in section [IT] we introduce the SVD method and test it against known input models in 
section IIHI We also show current constraints from the joint estimation of cosmological 
parameters and P{k) reconstructions in section HVl We discuss our results and conclude in 
section [V] 
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II. DIRECT INVERSION OF P k BY SINGULAR VALUE DECOMPOSITION 



Direct primordial power spectrum reconstruction requires the inversion of the following 
relations 



oo 

Cf Y = I ^-A?(k)Aj(k)P(k), (2) 
o 

where X and Y represent T, E, or B-type anisotropies, Cf Y are the angular power spectra 
for the XY combination and the Af(k) are the photon perturbation transfer functions. 
The transfer functions are obtained by integrating the full Einstein-Boltzmann system of 
differential equations Q, l^] . These describe the evolution of perturbations in the photon 
distribution functions in the presence of gravity and other sources of stress-energy. The 
functions determine all of the structure in the anisotropy spectra which arise after the 
initial conditions are set. Most notably the Cf contain distinct peaks due to the acoustic 
oscillation of the tightly coupled photon-baryon fluid in gravitational potential wells at the 
time of last scattering. The aim of any inversion method is to distinguish such features from 
any structure in the initial perturbation spectrum. 

For a finite sampling of the wavenumber space k Eq. (T5]) can be recast as an operator 
acting on the primordial spectrum P k 

Ct = '^^FgkPk, (3) 

k 

with operator 

F* Y = AlnkA? k Aj k , (4) 

where A In A; are the logarithmic k intervals for the discrete sampling chosen in the integration 
of the system of equations. 

A solution for P k cannot be obtained from a direct inversion of the F^ Y as it is numerically 
singular. This is due to the high level of degeneracy in the transfer functions relating the 
power at any wavenumber k to angular multipoles I. We instead approximate the inversion 
by using the SVD method, first reducing the degeneracy of the system and then inverting 
using the remaining orthogonal modes. 

The transfer functions can be factorised as 

Fik = 2j U u ,h ek ,V k , k , (5) 

i'k' 
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where the matrices U and V are unitary and of dimensions ri£ and n& respectively, and 
A is a non-negative, diagonal matrix with diagonal elements A&. For this application, the 
dimensions of the matrices are that ri£ < n& i.e. there are more equations than modes of 
interest. This results in some of the diagonal elements of A being singular (numerically zero) 
which prevent the inversion of the transfer matrix. 

The SVD method allows one to invert such a system by nulling the singular modes. This 
is achieved by creating an inverse A -1 where the diagonal elements are 1/Xk except where 
the value of is singular in which case it is replaced by 0. 

In practice we rank order the factorised modes in descending order of A& and all modes 
with condition number less than a threshold emax(Afc) are nulled, e ~ 0.038 for this work. 
Thus the method is a l k-to-€ compression of the system where we keep the least degenerate 
modes connecting the 3d Fourier space to the 2d angular multipole space. This is not to be 
confused with a signal-to-noise compression of the data which aims to select with respect to 



orthogonal modes of the covariance of the observations [481 ]. 

It is instructive to look at the first few orthogonal modes given by the columns of the 
U matrix. We plot the first six in Fig. [TJ These £-space modes are the least degenerate 
(or best determined) in the mapping provided by the CMB physics. In other words, in the 
absence of sample and noise variance, these modes pick out the Grange where observing 
the CMB will have the highest impact upon the reconstructed P(k). Not surprisingly, the 
first few modes are peaked around the angular scales where the acoustic signal from each 
polarisation combination is maximised. As these are the best constrained vectors in t space 
in the absence of all errors, they are not necessarily the basis vectors of P k which are most 
accurately constrained, [4^] show what they are for WMAP. We assume that this ordering 
of singular values is the optimal method for sorting the columns in U and V. 

Once A -1 has been computed the primordial power spectrum can be reconstructed by 
inverting a set of observed Ce using 



Pk — V kk iA k ,\,Uj, e C e , 

I'k'i 

« ( 6 ) 



Our choice of e ~ 0.038 is a conservative one with approximately 200 non-singular modes for 
a typical transfer matrix. Reconstructing the primordial power spectrum in this way means 
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FIG. 1: We plot the vectors in the matrix U associated with the 6 highest singular values. The 
vector corresponding to the highest singular value is shown in the top left panel, the second highest 
is shown in the top right panel, etc... These vectors are the modes best constrained in Ci in the 
absence of all sources of error. We decomposed F^i for a cosmology of f^/i 2 = 0.0226, f2 c /i 2 = 0.108, 
9 = 1.041 and r = 0.076. The blue/solid line is the TT mode, red/dashed is the TE mode and 
green/dotted is the EE mode. 
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that we are not using any degenerate modes which carry no information in k space but can 
increase the scatter in the reconstructed spectrum. However we are still susceptible to the 
scatter in the observed Cg since we have not used any noise weighting in this scheme. 

In practice we start our inversion process with a guess input spectrum, parameterised by 
the usual form Agk 71 " -1 . From this and our fiducial cosmological model we obtain a C™ odel 
spectrum. We use this to calculate the residual spectrum, 

£fres ^jobs £<model ^j^j 

so as to minimise the error induced in fc-space by the cut-offs in £, both on large and small 
scales. 

To remove high frequency oscillations in the data we apply a low-pass filter to the resultant 
Pfc. The following algorithm was used, 

plow-pass = ^ a) p k _^ (g) 

where a was taken to be 0.05. This method of smoothing leaves the first few points in the 
series strongly influenced by P\. Therefore one should take any effects seen at low k with 
a pinch of salt as these points are highly correlated. The covariance matrix was altered 
appropriately by a lower-triangular matrix representing this filter. 

We then proceed to bin the reconstructed power spectrum using the optimal binning 



method of 49j. This binning method estimates a series of ranges in k over which the signal- 
to- noise in the measured primordial power spectrum is constant. Many of the data points are 
highly correlated with their nearest neighbours and optimal binning gives a clear indication 
of the scales on which we have independent information. The centre of the k bins chosen by 
the cosmological tool CAMB [47j change when the input cosmological parameters are altered. 
This is a problem when it is run over many realisations as in the case of a Markov Chain. 
We choose the optimal binning method to find the standardised output of centres and sizes 
of k bins for each call of the CAMB routine. 

To find the optimal binning of the reconstructed power spectrum we investigate how the 
uncertainty in the Gt transfers into uncertainty in the primordial power spectrum. For this 
purpose, we need to define the primordial power spectrum as a series of top-hat bins: 

P(k) = Y,^B(k)Q B , (9) 

B 



where Qb is the power in each bin B and wb — 1 if k G B and otherwise. To obtain the 
errors for these bins we define the Fisher matrix for the Cg by 

Ma- = (SCw)' 1 , (10) 

where 5Cgg/ is the diagonal matrix of the squares of the variances in each measurement of 
Cg. To transfer the given errors from the Cg to other parameters we use the Jacobian, 

where, 6 a and 6p represent the bins of the primordial power spectrum. The derivative of the 
Cg with respect to the primordial power spectrum is the average radiation transfer function 
in each bin: 

8Cp f k max dk 



8P(*> J k B T ™'<*>- < 12 ' 
To calculate the signal-to-noise ratio in each bin we take the inverse square root of the 
diagonal elements of M af3 to be the noise and the amplitude of the primordial power spectrum 
to be the signal. We then arrange the bins to have the same signal-to-noise over our k range. 
Our algorithm will result in more bins where the signal-to-noise ratio is greater, sampling 
more finely where the signal is strongest. 

We construct a signal vector, S, which contains the amplitude of the primordial power 
spectrum for all the bins and weight our Fisher matrix by this vector 

- =S a M afi S p , (13) 

iV / a/3 

where there is no Einstein summation. The square root of the diagonal elements of this 
matrix are the S/N of the bins. 

We start our algorithm with the maximum number of bins possible in our k range. This 
is set by the usual properties of the Fourier transform. These imply that the scale of the 
survey not only determines k min , but also gives a lower bound upon the resolution, Ak min : 
narrower bins would become highly correlated. Therefore, we set up a series of bins with 
the properties 

£ min 2 A£ 1 

k min = —J— = ~T and {^k) min = — = — , (14) 

a a ua «A a a 

where cLa = 14.12Gpc (value given by WMAP5) is the angular diameter distance to the 

surface of last scaterring. We set k max = 0.08Mpc _1 as the reconstruction process is limited 



past this regime due to the cut-off in i at 1000. A Fisher matrix is then constructed for this 
set of bins and weighted by the signal vector. We choose a signal-to-noise value, that cannot 
be less than the maximum value seen in any bin, as our target ratio. The binwidths are 
increased in order to reach the target value at each bin. To obtain this 'optimal' binning we 
iterate until the bin with the smallest signal-to-noise ratio is within 5% of the target ratio. 

A. Reconstruction with cosmological parameter fitting 

The reconstruction method described above is fast and can be carried out at each ran- 
dom sample of a MCMC exploration of the cosmological parameter space. Inserting the 
reconstruction as part of an MCMC exploration we can account for the variance induced in 
the primordial power spectrum due to the dependence of the radiation transfer function on 
the cosmology. 

We have modified the cosmomc [|3(3] package by introducing the reconstruction at each 
chain evaluation using the inverse of the transfer function computed for each combination 
of parameters acting on the 'observed' CMB angular power spectrum. The reconstructed 
spectrum is then itself used to compute the final Ce which are used to calculate the likelihood 
at the chain step. 

In principle the chains would probe the reduced set of parameters; the physical densities of 
baryons Qbh 2 , and of cold dark matter Q c h 2 , the angular diameter distance parameter 9, and 
optical depth parameter r. The primordial power spectrum parameters n s and A s become 
irrelevant and need not be probed since the power spectrum is being reconstructed directly. 
However, in practice, we do include power law spectral parameters which determined the 
shape the template model (jTj) and we marginalise over the spectral parameters in order to 
account for any sensitivity of the reconstruction to the assumed C™ odel . 

The immediate advantage of combining the reconstruction with an MCMC method is 
that we can then calculate the variance in the resulting spectrum due to the random nature 
of the transfer function. We do this by including the binned amplitudes for the reconstructed 
spectrum as 'derived' parameters when analysing the chains. The covariance of the chains 
is then mapped into a covariance for the binned power spectrum. 

We also need to account for the variance due to the errors in the observed CMB data. 
This is not accounted for in the MCMC chains since we always use the same observed Ce 



9 



data to reconstruct the spectrum. In principle this contribution to the variance and that 
from the transfer function are correlated, however this is difficult to quantify without includ- 
ing MCMC steps over realisations of the observations. We therefore make a conservative 
estimate of the final error in the reconstructed spectrum by adding the covariance matrix ob- 
tained from the MCMC chain and that obtained by rotating the error matrix of the observed 
Ce as 

SPbb' = 2_s ^b£^u'Fb'\' (15) 

w 

where F^l is the bin- averaged contribution from Fj^ . 



III. APPLICATION OF THE SVD INVERSION 



A. Tests on simulated CMB data 



We start by testing the reconstruction algortihm on a set of input spectra with known 
features as in 181 ] . Fig. [2] shows that the method accurately reconstructs test features on 
the input spectra. We used a fiducial cosmological model of fl b h 2 = 0.0226, Q c h 2 = 0.108, 
6 = 1.041 and r = 0.076. The starting spectrum used to obtain Cg model is shown by the 
black, long-dashed line. We assume there are no errors on the input C° bs (between 1 = 2 
and C = 1000) and observe that the reconstructed Pk generally picks out the input features. 
We limit our method to fit between k = 0.0013 and k = 0.08 as this is the range with the 
highest signal-to-noise in the WMAP data. The features chosen are the same as in [3] 
to allow comparisons with this method. They are a standard power law with n s = 0.96 
but with an amplitude 90% of that of the best-fit WMAP model, a running n s model with 
dn.Jdlnk = —0.037, a power law with a sharp, compensated feature at k = 0.02 Mpc -1 



16, 51 
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151 ] and a power law with superimposed sinusoidal oscillations 



All four features are clearly recovered to varying degrees when using TT, 



TE or EE, however we find a phase offset between the reconstructed ad input spectrum 
for the case where the input spectrum includes oscillations as in the lower panels of Fig. [21 
The offset is stable with respect to the presence of the smoothing kernel, the number of 
singular values cut from the inversion, and with repsect to the number of k bins and range. 
The reconstructed Pk given by both TE and EE contains 'glitches' not present in the TT 
reconstruction. These regions correspond to regions where there is little information int he 
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FIG. 2: The reconstruction of several test spectra. The test models used to generate the simulated 
Ci, shown in black (thick/solid) curves, are (a) A 10% decrease in power from the WMAP5 best fit 
amplitude, (b) the WMAP5 best fit model including running dn s /dlnk = —0.037, (c) a localised 
feature at around k = 0.02 Mpc -1 , and (d) a model with sinusoidal oscillations superimposed 
on the best fit power law spectrum. The black (long-dashed) curves show the best fit spectrum 
used to minimise any cut-off effects at the ends of the I regions. The blue (solid) curves are 
the reconstructions using total intensity data whereas the red (dashed) curves and magenta (dot- 
dashed) curves use TE and EE data respectively. The Cg forecasts assumed an experiment with 
no noise and an ^ max = 1000. 
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FIG. 3: The reconstruction of several spectra when the incorrect parameters are used for the 
fiducial cosmological model. The black (solid) line is the correct input spectrum. In each panel 
we change a single parameter in the the fiducial model (f^/i 2 = 0.0226, Q c h 2 = 0.108, 8 = 1.041 
and r = 0.076). The red (dot-dashed) line shows the effect of changing a parameter by la in a 
positive direction from the WMAP best fit model, the magenta (dashed) line shows the effect of a 
la shift in a negative direction. The top left panel shows the effect of varying f^/i 2 , top right of 
Q c h 2 , bottom left of 6 and the bottom right of r. 

TE and EE spectra and the reconstruction is still degenerate. 

The cosmological parameters obtained from a traditional Markov Chain Monte Carlo 
search are used to give us our fiducial operator, F^. The parameters however have errors 
present upon them, which are not usually incorporated into errors on the final reconstructed 
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FIG. 4: We show the unbinned errors on a reconstructed P(k) for a forecasted Planck dataset. The 
cosmological parameters have been fixed. The blue (solid) lines show the la confidence regions 
obtained from TT measurements, with red (dashed) and magenta (dot-dashed) showing the same 
lcr bounds for TE and EE respectively. 

Pfc. If inaccurate parameters have been used to calculate the fiducial operator very spe- 
cific signatures would be expected to show up in the reconstructed from of P^. We show 
in Fig. [3] how these signatures appear in the reconstructed for TT anisotropics for four 
cosmological parameters, Q^h 2 , Q c h 2 , h and r. If any features with this form are observered 
one should attribute this to an incorrect estimation of the parameters and not some funda- 
mental physics. A check for this would be if the corresponding features also show up in the 
reconstructed TE and EE spectra. 
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FIG. 5: Predicted constraints on Q^h 2 , Q c h 2 , 6 and r from Planck, when is given total freedom. 
The solid red vertical lines indicates the input values for each of the parameters. The black solid 
curves show the marginalised probability distribution. 

We have also tested the inversion on simulated CMB data with similar experimental 
properties as the recently launched Planck satelite mission [53]. We assume a total of 12 
detectors with NETs of 64/aK/ ^/s observing 80% of the sky over 12 months with a resolution 
of 7 arcminutes FWHM. We calculate errors around our fiducial CMB best-fit models in both 
total intensity and polarisation spectra for this experimental setup and use these together 
with Ci samples on the fiducial model to test the inversion method's properties. We consider 
multipoles of I < 1000 for both total intensity spectra and polarisation. We have not taken 
into account any residual error from foreground subtraction in our forecasts. Thus our 
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FIG. 6: Predicted constraints on the form of from Planck. The red line is the reconstructed P^ 
at the best-fit point. We show in black the marginalised values of each bin, the error bars represent 
the lo" error. 

forecast are on the optimistic side of the accuracy achievable in the case of polarisation where 
foreground removal will certainly have a significant impact on errors at i < 1000. In the 
case of total intensity spectra we are significantly underestimating the accuracy achievable 
by Planck as we expect to obtain well measured Cgs well above £ of 1000. We do not 
consider these modes as it significantly increases the time required to perform the SVD. 
Consideration of total-intensity modes past i of 1000 will increase the accuracy obtainable 
for Planck on it also allows one to probe k past 0.08. It would be desirable to perform 
this process with a greater I range when the Planck data becomes available. 
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To compare the accuracy obtainable with each of the anisotropy types we show the degree 
to which they each reconstruct a simple power law Pk in Fig. HI We use the same fiducial 
cosmological model as in Fig. [2] with a standard power law input Pk equal to the spectra 
being reconstructed. All the C° bs were placed on the fiducial model. Over the whole range 
of k total-intensity modes best constrain P^. But it is also true that there are regions in 
/c-space, for example between k = 0.022 and k = 0.035, where considering only the TE 
measurements can give us a more accurate estimation of Pk- EE measurements approach 
the accuracy of TT at a number of points, there are however regions where the errors become 
so large that any reconstructed Pk is meaningless (these correspond to the troughs of the 
EE spectrum). The regions of £-space corresponding to the most accurately measured k 
regions are the peaks of the spectra. Both TE and EE spectra are not significantly affected 
by changes to Pk at very low k, so errors in this region are artificially small. 

We also tested the combined P(k) and parameter estimation MCMC search as described 
in the previous section. The optimal binning method found 128 bins for the Planck exper- 
iment, where our target signal-to-noise value in each bin is 10. We choose this value so as 
to have approximately the same number of bins across the range 0.01 < k < 0.03 where 
WMAP best probes P k . 

In Fig. [5] we show the results of a cosmological parameter estimation from this process 
for a simulated Planck experiment. We find that the input parameters were accurately 
recovered after this process. We show the final reconstructed Pk in Fig. [6l The red line 
is the reconstructed Pk at its best-fit point. The error bars we show are obtained from 
combining the error from the marginalised distributions with the reconstruction errors given 
the observed CMB data. The errors are centred around the mean of the marginalised 
distribution for each bin. It is important to note that the errors are highly correlated. 
This explains the reduced scatter in the mean values compared to the size of the plotted 
errors. As was seen in Fig. [3] changing a single parameter by a small amount (in the manner 
a MCMC search does) creates a very specific signature on Pk for each parameter, where 
changing it slightly has a correlated effect upon the whole range of k. This explains the very 
high correlations observed accross the whole k range. It is at odds with the errors on P(k) 
associated with those of Gg_ at any best-fit point which are not correlated across large ranges 
of k. 
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FIG. 7: The current constraints on f^/i 2 , Q c h 2 , 9 and r, when P^ is given total freedom. The 
solid red vertical lines indicates the WMAP only best fit values when P(k) is parameterised in 
the usual fashion, the black solid curves represent the marginalised probability distribution of the 
WMAP only data. The blue dotted line shows the marginalised probability distribution of WMAP 
including the other datasets. 

IV. CONSTRAINTS FROM CURRENT CMB OBSERVATIONS 



We used two sets of currently available data to estimate the cosmological parameters in 
conjunction with a free unparameterised The first is the WMAP 5-year data alone 541 ] . 
in the second we combine this with that of SNIa, HST and BBN 55j, |56|, 157 



58, E]. 
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FIG. 8: The current constraints on the form of from WMAP data only. The red lines is the 
reconstructed at the best-fit point. We show in black the marginalised values of each bin, the 
error bars represent the la error. 

second set combination chosen because the non-CMB sets do not depend upon the form of 
Pk and can therefore give us independent and tighter constraints on most of the cosmological 
parameters. We run this data through the MCMC tool C0SM0MC in the same manner as we 
did for the test Planck data. 

In Fig. [7] we show the current constraints on f^/i 2 , Q c h 2 , 9 and r. The solid red vertical 
lines indicates the WMAP only best fit values when P(k) is parameterised in the usual 
fashion, the black solid curves represent the marginalised probability distribution of the 
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FIG. 9: The current constraints on the form of P^ from WMAP and the other datasets. The red 
lines is the reconstructed P& at the best-fit point. We show in black the marginalised values of 
each bin, the error bars represent the la error. 

WMAP only data. There is no disagreement between this and the WMAP bestfit model. 
The blue dotted line shows the marginalised probability distribution of WMAP including 
the other datasets. Here we observe some tension at around the la level in Q c h 2 . The 
inclusion of this data does not move the position of the likelihood peak significantly for each 
parameter. 

The optimal binning method gave us 86 bins in our k range for WMAP. We used the 
minimum Ak as described previously, this gave a target signal-to-noise ratio of 6. The 



19 



binned reconstructed Pk are shown in Figs. [8] and [9] for each of the two data sets. As in 
the Planck case, the results of the MCMC are highly correlated across the whole k range, 
in contrast to a single reconstruction which is not. We have similarly added the errors from 
the MCMC to a single reconstruction in each bin in quadrature. The red lines are the 
reconstructed Pk at their best-fit points. They error bars are centred around the best-fit 
marginalised value of each bin. For the WMAP only run we find that the limiting errors on 
Pk in the range 0.0075 < k < 0.05 come from uncertainty in the cosmological parameters, 
whereas this limiting range is around 0.0075 < k < 0.04 when we include the other datasets. 
No significant deviation from the standard power-law case is observed in either run. 

V. DISCUSSION 

The SVD based reconstruction method we have outlined provides is fast enough to be 
incorporated into full MCMC parameter fitting runs. We have tested the method and 
shown that it recovers the overall features of input spectra. We have applied the method to 
forecasted Planck data and current WMAP 5-year results. These results have allowed us to 
consistently combine the reconstruction with a full exploration of the parameter likelihoods 
for the first time. We have seen that the limiting factor in constraining the primordial 
spectrum over a large range of wavenumbers k comes from the uncertainty in cosmological 
parameters. Any claims of a detection of a feature in Pk must necessarily confront the 
degeneracy with the cosmological parameter space . This effect will be less important when 
the Planck data is released, however it must still be considered as the unprecedented accuracy 
offered by future data may lead to premature claims of a detection of an interesting feature. 

We observed some tension between the WMAP only best-fit model, when P(k) is param- 
eterised with the normal amplitude and tilt, and our method in the marginalised probability 
distribution of Q c h 2 . They disagree at around the la level which is not overly significant, 
however it is an indication that there may be some departure from a standard power law in 
P{k). Planck will certainly determine if this is the case. 

There are other currently available CMB datasets, which could expand the range of k 
probed. However expanding the method to include multiple data sets is a non-trivial task 
due to the binning of the Ce data for sub-orbital experiments and the increased running 
time required for the expanded fc-range. We leave this analysis for future work. 
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In the future, as CMB polarisation data becomes increasingly accurate, it will be desirable 
to perform a joint inversion of total intensity data along with polarisation data. It is not 
clear how to extend the SVD based method to include all polarisation modes simultaneously 
since a HO SVD (Higher-Order SVD) step would probably be required. On the other hand 
this would give the best estimate of P^ given any dataset and would help to reduce the 
correlations found in the reconstructed Pk by increasing the degrees of freedom that can be 
effectively constrained. 
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